{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# General import \n",
    "\n",
    "import numpy as np\n",
    "import scipy.sparse as sparse\n",
    "from scipy.integrate import ode\n",
    "import time\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pyMPC import\n",
    "\n",
    "from pyMPC.mpc import MPCController"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System dynamics ##\n",
    "\n",
    "The system to be controlled is an inverted pendulum on a cart (see next Figure).  \n",
    "\n",
    "<img src=\"img/cart_pole.png\" width=\"250\" align=\"center\"/>\n",
    "\n",
    "The system is governed by the following differential equations:\n",
    "\n",
    "\\begin{equation}\n",
    " \\begin{aligned}\n",
    " (M+m)\\ddot p + ml\\ddot\\phi \\cos\\phi - ml \\dot \\phi ^2 \\sin \\phi + b\\dot p &= F \\\\\n",
    " l \\ddot \\phi + \\ddot p \\cos \\phi - g \\sin\\phi &= -f_\\phi\\dot \\phi\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "Introducing the state vector $x=[p\\; \\dot p\\; \\phi\\; \\dot \\phi]$ and the input $u=F$, the system dynamics are described in state-space by a set of an nonlinear ordinary differential equations: $\\dot x = f(x,u)$ with\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "  f(x,u) &= \n",
    " \\begin{bmatrix}\n",
    " x_2\\\\\n",
    " \\frac{-mg \\sin x_3\\cos x_3 + mlx_4^3\\sin x_3 + f_\\phi m x_4 \\cos x_3 - bx_2 + u }{M+(1-\\cos^2 x_3)m}\\\\\n",
    " x_3\\\\\n",
    " \\frac{(M+m)(g \\sin x_3 - f_\\phi x_4) - (lm x_4^2 \\sin x_3 - bx_2 + u)\\cos x_3}{l(M+(1-\\cos^2 x_3)m)}\n",
    " \\end{bmatrix}\\\\ \n",
    " \\end{split}\n",
    " \\end{equation}\n",
    "\n",
    "For MPC control design, the system is linearized about the upright (unstable) equilibrium point, i.e., about the point $x_{eq} = [0, \\; 0\\;, 0,\\; 0]^\\top$.\n",
    "The linearized system has form $\\dot  x =  A_c x + B_c u$ with\n",
    "\n",
    "\\begin{equation}\n",
    "  A = \n",
    " \\begin{bmatrix}\n",
    " 0& 1& 0& 0\\\\\n",
    " 0& -\\frac{b}{M}& -g\\frac{m}{M}& f_\\theta\\frac{m}{M}\\\\\n",
    " 0&0&0&1\\\\\n",
    " 0&\\frac{b}{Ml}& \\frac{g(M+m)}{Ml}&-\\frac{(M+m)f_\\theta}{M l}\n",
    " \\end{bmatrix},\\qquad B= \n",
    " \\begin{bmatrix}\n",
    " 0\\\\\n",
    " \\frac{1}{M}\\\\\n",
    " 0\\\\\n",
    " -\\frac{1}{Ml}&\n",
    " \\end{bmatrix}\n",
    " \\end{equation}\n",
    " \n",
    "Next, the system is discretized with sampling time $T_s = 50\\;\\text{ms}$. Here we just use a Forward Euler dsicretization scheme for the sake of simplicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constants #\n",
    "\n",
    "M = 0.5\n",
    "m = 0.2\n",
    "b = 0.1\n",
    "ftheta = 0.1\n",
    "l = 0.3\n",
    "g = 9.81\n",
    "\n",
    "Ts = 50e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# System dynamics: \\dot x = f_ODE(t,x,u)\n",
    "\n",
    "def f_ODE(t,x,u):\n",
    "    F = u\n",
    "    v = x[1]\n",
    "    theta = x[2]\n",
    "    omega = x[3]\n",
    "    der = np.zeros(4)\n",
    "    der[0] = v\n",
    "    der[1] = (m * l * np.sin(theta) * omega ** 2 - m * g * np.sin(theta) * np.cos(theta) + m * ftheta * np.cos(theta) * omega + F - b * v) / (M + m * (1 - np.cos(theta) ** 2))\n",
    "    der[2] = omega\n",
    "    der[3] = ((M + m) * (g * np.sin(theta) - ftheta * omega) - m * l * omega ** 2 * np.sin(theta) * np.cos(theta) - (F - b * v) * np.cos(theta)) / (l * (M + m * (1 - np.cos(theta) ** 2)))\n",
    "    return der"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Linearized System Matrices \n",
    "    \n",
    "# Continuous-time system matrices, linearized about the upright, unstable equilibrium point \n",
    "Ac =np.array([[0,       1,          0,                  0],\n",
    "              [0,       -b/M,       -(g*m)/M,           (ftheta*m)/M],\n",
    "              [0,       0,          0,                  1],\n",
    "              [0,       b/(M*l),    (M*g + g*m)/(M*l),  -(M*ftheta + ftheta*m)/(M*l)]])\n",
    "\n",
    "Bc = np.array([\n",
    "    [0.0],\n",
    "    [1.0/M],\n",
    "    [0.0],\n",
    "    [-1/(M*l)]\n",
    "])\n",
    "\n",
    "[nx, nu] = Bc.shape # number of states and number or inputs\n",
    "\n",
    "# Simple forward euler discretization\n",
    "Ad = np.eye(nx) + Ac*Ts\n",
    "Bd = Bc*Ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MPC reference input and states (set-points)\n",
    "\n",
    "xref = np.array([0.3, 0.0, 0.0, 0.0]) # reference state\n",
    "uref = np.array([0.0])    # reference input\n",
    "uminus1 = np.array([0.0])     # input at time step negative one - used to penalize the first delta u at time instant 0. Could be the same as uref.\n",
    "\n",
    "# Constraints\n",
    "xmin = np.array([-100.0, -100, -100, -100])\n",
    "xmax = np.array([100.0,   100.0, 100, 100])\n",
    "\n",
    "umin = np.array([-20])\n",
    "umax = np.array([20])\n",
    "\n",
    "Dumin = np.array([-5])\n",
    "Dumax = np.array([5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MPC constraints\n",
    "xmin = np.array([-100.0, -100, -100, -100])\n",
    "xmax = np.array([100.0,   100.0, 100, 100])\n",
    "\n",
    "umin = np.array([-20])\n",
    "umax = np.array([20])\n",
    "\n",
    "Dumin = np.array([-5])\n",
    "Dumax = np.array([5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MPC cost function weights\n",
    "\n",
    "Qx = sparse.diags([0.3, 0, 1.0, 0])   # Quadratic cost for states x0, x1, ..., x_N-1\n",
    "QxN = sparse.diags([0.3, 0, 1.0, 0])  # Quadratic cost for xN\n",
    "Qu = 0.0 * sparse.eye(1)        # Quadratic cost for u0, u1, ...., u_N-1\n",
    "QDu = 0.01 * sparse.eye(1)       # Quadratic cost for Du0, Du1, ...., Du_N-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<scipy.integrate._ode.ode at 0x2024e539b70>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Initialize simulation system\n",
    "\n",
    "phi0 = 15*2*np.pi/360\n",
    "x0 = np.array([0, 0, phi0, 0]) # initial state\n",
    "t0 = 0\n",
    "system_dyn = ode(f_ODE).set_integrator('vode', method='bdf')\n",
    "system_dyn.set_initial_value(x0, t0)\n",
    "system_dyn.set_f_params(0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "Np = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize and setup MPC controller\n",
    "\n",
    "K = MPCController(Ad,Bd,Np=Np, x0=x0,xref=xref,uminus1=uminus1,\n",
    "                  Qx=Qx, QxN=QxN, Qu=Qu,QDu=QDu,\n",
    "                  xmin=xmin,xmax=xmax,umin=umin,umax=umax,Dumin=Dumin,Dumax=Dumax)\n",
    "K.setup() # this initializes the QP problem for the first step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate in closed loop. Use MPC model as real system\n",
    "\n",
    "# Simulate in closed loop\n",
    "[nx, nu] = Bd.shape # number of states and number or inputs\n",
    "len_sim = 10 # simulation length (s)\n",
    "nsim = int(len_sim/Ts) # simulation length(timesteps)\n",
    "xsim = np.zeros((nsim,nx))\n",
    "usim = np.zeros((nsim,nu))\n",
    "tsim = np.arange(0,nsim)*Ts\n",
    "\n",
    "time_start = time.time()\n",
    "\n",
    "t_step = t0\n",
    "uMPC = uminus1\n",
    "for i in range(nsim):\n",
    "    xsim[i,:] = system_dyn.y\n",
    "\n",
    "    # MPC update and step. Could be in just one function call\n",
    "    K.update(system_dyn.y, uMPC) # update with measurement\n",
    "    uMPC = K.output() # MPC step (u_k value)\n",
    "    usim[i,:] = uMPC\n",
    "\n",
    "    # System simulation\n",
    "    system_dyn.set_f_params(uMPC) # set current input value\n",
    "    system_dyn.integrate(t_step + Ts)\n",
    "\n",
    "    # Time update\n",
    "    t_step += Ts\n",
    "time_sim = time.time() - time_start\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAJOCAYAAAB1IEnpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl8VNX9//HXZ7KSBAgmEJYEwiYQlL1iRSvUpWoVW78o2Lrw0xbF7Wtta7XWLrS2tmrdKm61X1Gr1Fpb0bpVa1wqpQKCEgiWsIQICAkEspD9/P6YScwGhOSGO8m8nw/nkZl77vJJzpC8PffOueacQ0REREQ6LuB3ASIiIiLdhYKViIiIiEcUrEREREQ8omAlIiIi4hEFKxERERGPKFiJiIiIeETBSkTCjpkNNrNSM4s6yDqlZjask47/KzO7vp3bXmdmt3tdk4h0DaZ5rESko8xsM5AG1AJlwMvAtc65Uo/2nw085Zz7vRf7O8Sx+gKrgBHOuf3t2D4e2ABMcs7t9Lo+EQlvGrESEa+c45xLAiYBXwB+5HM97TUXeLk9oQrAOVcBvAJc4mVRItI1KFiJiKecc58SDBbHAJjZQDNbYma7zWyDmX27fl0zO87MlpvZPjP7zMx+G1qeaWbOzKLN7DbgJOB3odN/vwut48xsROh5bzN7wsx2mdkWM/uRmQVCbXPN7D0zu9PM9pjZJjM78yDfwpnA241qnG5mBWZ2o5ntNLPtZvY1MzvLzD4JfV8/bLaPbOCrHfxRikgXFO13ASLSvZhZBnAW8Hxo0TNADjAQGA38w8w2OufeBO4F7nXOPWlmSYTCWGPOuVvMbBoHPxV4P9AbGAakAK8D24HHQu1TgUVAKjAPeMzMBrnWr4U4FljfbFl/IB4YRHBE61HgH8BkYDCwwswWO+c2htZfB4w/QK0i0o1pxEpEvPI3MysG3iM44vPLUMg6EfiBc67CObcK+D1wcWibamCEmaU650qdc/8+3IOGLnCfDdzsnCtxzm0G7mp0DIAtzrlHnXO1BAPWAILXhLUmGShptqwauM05Vw0sJhjQ7g0dL4dgcBzXaP0SgkFPRCKMgpWIeOVrzrlk59wQ59xVoWuUBgK7nXONg8oWgiM/AJcDRwO5ZvaBmZ3djuOmArGh/bZ2DIAd9U+cc+Whp0kH2N8eoGezZUWhUAZQf+3VZ43a9zfbX09g7yErF5FuR8FKRDrTNuAoM2scVAYDnwI45/7rnLsQ6Af8GnjOzBJb2c/BPr5cSHBEaUhrx2iHjwiGvY4YA6zu4D5EpAtSsBKRTuOc2wq8D/zKzOLNbBzBUao/ApjZRWbW1zlXBxSHNqttZVefEbx+qrVj1ALPAreZWU8zGwLcADzVzrJfBk5u57b1TiZ4Ab+IRBgFKxHpbBcCmQRHr/4K/MQ5949Q2xlAjpmVEryQfU5ouoLm7gVmhT7Vd18r7dcSnD9rI8FrvJ4G/tDOep8AzjKzHu3ZODSP1VkEr+USkQijCUJFRJoxs18CO51z97Rj22uBDOfcjd5XJiLhTsFKRERExCM6FSgiIiLiEQUrEREREY8oWImIiIh4xLdb2qSmprrMzMxOPUZZWRmJia1NiSN+Ur+EJ/VLeFK/hCf1S3jqzH5ZsWJFoXOu76HW8y1YZWZmsnz58k49RnZ2NtOnT+/UY8jhU7+EJ/VLeFK/hCf1S3jqzH4xsy2HXkunAkVEREQ8o2AlIiIi4hEFKxERERGPKFiJiIiIeMS3i9fFW8459uzZw44dOygpKaGsrIyamhri4uKaPJKSkujfvz+xsbF+lywiItLtKFh1QXv27GHp0qUsXbqUnJwc1q1bx6ZNm6isrGzzPvr27cugQYMYOHAg6enpjBo1ilGjRjF69GgyMzOJiorqxO9ARESke1Kw6iLWr1/PX/7yF55//nlWrlyJc46oqChGjBhBVlYWZ599NgMHDqR///707t2bxMREoqOjqaysbHhUVFRQUlLCtm3bGh6ffvopy5Yto6ioqOFYsbGxjBw5kvHjxzNx4kQmTJjAhAkTSE1N9fEnICIiEv4UrMJYdXU1f/vb37jvvvt47733AJg6dSo//elPOemkk5g6dSoJCQmeHKuoqIj169eTm5tLbm4u69at45133uHpp59uWGfQoEENIWvixIlMnDiRoUOHYmae1CAiItLVKViFIecczz77LDfffDObNm1i2LBh3HHHHcyZM4f09PROOWZKSgonnHACJ5xwQpPlhYWFrF69mlWrVjU8Xn31VWprawFITk5uCFmTJk1i4sSJjBo1SqcSRUQkIilYhZnVq1dzxRVXsGzZMsaPH88LL7zAV7/6Vd+CSmpqKqeccgqnnHJKw7KKigrWrFnDypUr+fDDD1m5ciULFy6koqICgISEBMaNG9cQtCZNmsTYsWOJi4vz5XsQERE5UhSswkRdXR333nsvN910E3369OEPf/gDl1xySViO/MTHxzNlyhSmTJnSsKympobc3NwmYeupp55i4cKFAMTExDB27FgmTZpEUlISsbGxjB8/XvfaEhGRbkXBKgzs27ePOXPm8MorrzBz5kwee+yxLneheHR0NMcccwzHHHMMl1xyCRAMixs3bmwIWitXrmTJkiUUFhZy3333YWYNF9+PHTuWrKwssrKyGD16ND169PD5OxIRETl8ClY+27FjB2eeeSZr1qzhgQceYP78+d3mYvBAIMCIESMYMWIE559/PhC8fuy5554jPj6elStXsmbNGnJycvj73/9OTU0NAGbGsGHDyMrKYuTIkQ37GD58OIMHDyY6Wm9bEREJT/oL5aO8vDxOO+00PvvsM1588UXOOOMMv0vqdGZG3759mT59Ouecc07D8qqqKjZs2MDatWvJyclh7dq1rF27ljfeeIP9+/c3rBcdHc3QoUMZPnw4w4YNIyMjo8lj4MCBupZLRER806ZgZWZnAPcCUcDvnXO3N2u/ErgaqAVKgXnOubUe19qtbNu2jVNPPZWSkhLeeustjjvuOL9L8lVsbGzDqcBZs2Y1LK+rq2P79u3k5eWxYcOGhkdeXh7//ve/KS4ubrGvtLQ0MjIyGDBgAGlpafTr149+/fo1PK//mpycTExMzJH8NkVEpJs7ZLAysyjgAeA0oAD4wMyWNAtOTzvnHgqtPxP4LdD9h1/aqbi4mDPOOIPCwkKys7OZPHmy3yWFrUAgwKBBgxg0aBBf+tKXWrSXlpZSUFDA1q1b2bp1a5Pn+fn5LF++nJ07dzZMD9FcYmIiycnJB30kJSWRkJBwyEdcXBwxMTEEAroFp4hIpGrLiNVxwAbn3EYAM1sMnAs0BCvn3L5G6ycCzssi22vC9ddDcnLThRdcAFddBeXlcNZZLTeaOzf4KCyERiMnDebPh9mzYetWuPjilu3f/S6ccw6sXw9XXNGiuerGG5l5++3ErVtHflYWfb773aYr/PKXcMIJ8P778MMfttz/PffAhAnwxhvwi1+0bH/4YRg1Cl58Ee66q2X7k09CRgb86U/w4IMt2597DlJT4fHHg4/mXn4ZEhJg4UJ49tmW7dnZwa933gkvvdS0rUcPeOWV4POf/xzefLNpe0oK/OUvwec33wxLlzZtT0+Hp54KPr/+eli1iiRgdOjB0UfDI48E2+fNg4oKSE7GHX00NdXV7Bs2jI8vv5zPPvuM8XfeSXxhITU1NcHHZ5/xUWkpd+3fT25uLndv3Uqvqqomb+Q3gfqf+MtAebNv/SXgLiAqKop/1tVhgQBmRsAMCwR4JSmJZ1NT6RUdzcItWwiYgRkGYMaraWm8PnAgvaur+VlODoSutatvf2XIEN4bNIh+lZV8Z+XKz7cNWTJyJMsHDGDAvn1cuWoVuM+rd8CzI0eyKjWVoXv38u2cnKb/SJ3jqf79uX3wYMbs2cPc//63YXm9B0aOZENSEpN27+aizZsb9lvvjuHD2ZqQwAlFRcz59NMWvwR+Nnw4O+PiOKWwkPM++6xF+03Dh1McHc3ZhYWcXVjYYv/XDR/OfjPOLyzk9D17WrRfPmwYAJfs2sXJ+z7/leSASjOuGDwYgCt37eL4srIm7cVRUVw7YAAA3y0sZEKj088A26Oj+W6/fjjnuLWoiDHNbh+1MTqam1NSAPhVURFDQ9cL1lsbG8uCPn0AuKewkP7Ngv7KuDh+E/pd9dCuXSTX1TW0xdbVcVdCAvf17g3Aop07iXNNf3r/7NGDR3r1AmDxZ5/R3N8TEniyZ0/i6+p4fNeuFu3PJSbyXFISfWpreTD0s2/sqaQkXkpMZEBNDXc3uktDvUd79uTNhASGVVfzy927W7Tf37s3/4qPJ6uqih+H+q6xO5KTWREXx+TKSr7fyij0gj59WBsby7SKCq7du7dF+w+POoqNMTGcUl7Ot0tKWrR/JyWF7dHRnF1WxkWlpS3a56emsicqilmlpcxq9N6oN7dvXyoCAS4uKeGr5cF/+bF1dSwN/U/UnLQ0AObt28eXm713Ks24tF8/AK7bu5cTQtPS1CsOBLiyb18AbiwuZlKz99aOqCiuD32Y6cd79pBVVdWkfVMnvvcA3o+P7xLvvZtvvpkrr7yyRbsf2hKsBgFbG70uAKY2X8nMrgZuAGKBL7e2IzObB8yD4Oma7Po/wp3k2NraFqeKdn7yCduyswlUVDCulX/AO3Jz2ZGdTczevYxtpf3TnBx2ZWcTt3MnY1pp3/rxxxT17EmP/HxGtdL+21/8gneXLuWBb38b+89/WtS3ceVK9lVV0WvNGoa1sv2G5cspLS6mz+rVDGmlff2yZezfvp2Ujz8mo5X2dUuXUpmXR9+cHAa10p7zr39R3bs3/XNz6d9K+0fvvENdfDwDP/mEfq20rwr1aUZeHinN2mv37+fj7GxKS0vZtGkTfZq1V9fVkRPafmh+Pr2btVfGxLAu1D6ioICkZu3l27bxSaj96G3bSGjWXn+tVlpaGn2Sk4mrrm7S3nPsWAZ++9sAjP3xj4net4+6urqGR69Roxh4+ulUVlZy9IMPElVZGWxzjrq6Oo5LT2duVhY1NTWkvvYazrkmj9TUVAYMGEB0VRUBM5xz4Bx1jerbvXs3tdXVVIV+eTpoCDfbd+xgfWkpxVVVlNb/8m/0Sy4vL4+V27czrLqa0kZ/POrDV35+PuuLighUVlIe+lk0DmalpaVs2bKFo/bvp7KykuYfodi9ezeFFRXsKy2luv6DBs1+vmXOUVlZSW3ol3Pj9traWmpqaqirq2sRqozg6GRUVBSBUCBt2D70PD4+HouKIjYmhkBUVIv6kkN/HBJKS4lu9sfNBQL0C/1xS9y/n9j6vg/tu0dMDBkZGQD0rKoizrkm+0+Kj2f48OEA9K6tpUejP94GJCcmkpWVBUCfNWtIbPbHOaVXL7LGjAluv3o1ic3+uKYmJ5M1alTw+B9+SGKjP57OOfqmppI1YkSw/vJy4pr9cUzr14+soUOD7a0Ei/4DBpA1eDBxtbUtagMYOHAgWenp9KqqIrG8+f8yQHp6OlkDBpC6fz+JzX62AIMHDyarXz8GlZa2+N4AhgwezJ7UVIbt20fiunUt2jMzM9nfpw+Ze/aQ+MknLdqHDRsGvXoxpLCQxLy8Fu3Dhw8nPimJwTt3krhpU4v2kSNH0qdHD9K3bycxP79F+6hRo9gXG8vAggISP/20RfuYMWOojIqif34+idu3A8F+qX+f1vd92qZNJO7c2WTb6Kiohva+GzaQ2Cwc1IYugwBIXb+exGa/t3rHxze0p6xbR+K+fU3a+3Tiew+gb0pKl3jvFRcXkx36+9LZ2eJQzLmDDy6Z2fnAV5xz3wq9vhg4zjl37QHW/0Zo/UsPtt8pU6a45cuXt6/qNsrOzmb69OmdeozDsWjRIubOncuNN97Ir3/9a7/L8U249YsEqV/Ck/olPKlfwlNn9ouZrXDOTTnUem25GKQAyGj0Oh3YdpD1FwNfa8N+I8qqVau48sormTFjBrfddpvf5YiIiEgnaEuw+gAYaWZDzSwWmAMsabyCmY1s9PKrwH+9K7HrKy8vZ86cOaSkpLB48WLNwyQiItJNHfIvvHOuxsyuAV4jON3CH5xzOWa2AFjunFsCXGNmpwLVwB7goKcBI833v/991q9fz5tvvtlwnYeIiIh0P20aOnHOvUzwg1CNl/240fP/9biubuPll19m4cKF3HDDDXz5y61e0y8iIiLdhCbc6URFRUVcdtllHHvssbquSkREJALoYp9OdMMNN1BUVMRrr71GfHy83+WIiIhIJ9OIVSd5/fXXeeKJJ7jpppsYP3683+WIiIjIEaBg1QnKysq44oorGDVqFLfccovf5YiIiMgRolOBneDWW29l8+bNvPPOOzoFKCIiEkE0YuWx//znP9x7773Mnz+fk046ye9yRERE5AhSsPJQdXU13/rWtxgwYAC/+tWv/C5HREREjjCdCvTQHXfcwccff8wLL7xA79DdwEVERCRyaMTKI5988gkLFizg/PPPZ+bMmX6XIyIiIj5QsPKAc47rrruOuLg47rvvPr/LEREREZ/oVKAHXnjhBV577TXuuece+vfv73c5IiIi4hONWHXQ/v37uf766znmmGO4+uqr/S5HREREfKQRqw769a9/zZYtW3jrrbeIjtaPU0REJJJpxKoDduzYwW9+8xsuuOACpk+f7nc5IiIi4jMFqw647bbbqKqq4rbbbvO7FBEREQkDClbttHnzZh5++GEuv/xyRowY4Xc5IiIiEgYUrNppwYIFBAIBbr31Vr9LERERkTChYNUO69evZ9GiRVx11VWkp6f7XY6IiIiECQWrdrj33nuJiYnhpptu8rsUERERCSMKVoeppKSEJ598ktmzZ9OvXz+/yxEREZEwomB1mJ5++mlKS0uZP3++36WIiIhImFGwOgzOOR588EHGjx/P1KlT/S5HREREwoymCj8My5YtY/Xq1Tz00EOYmd/liIiISJjRiNVhePDBB0lKSuIb3/iG36WIiIhIGFKwaqP9+/fz3HPP8Y1vfIOePXv6XY6IiIiEIQWrNsrOzqa8vJzzzjvP71JEREQkTClYtdFLL71EYmIiJ598st+liIiISJhSsGoD5xwvvfQSp512GvHx8X6XIyIiImFKwaoN1qxZQ35+Pl/96lf9LkVERETCmIJVG/z9738H4KyzzvK5EhEREQlnClZt8NJLLzF58mQGDhzodykiIiISxhSsDqGwsJClS5fqNKCIiIgckoLVIbz22mvU1dVx9tln+12KiIiIhDkFq0N45513SE5OZvLkyX6XIiIiImFOweoQli1bxnHHHUcgoB+ViIiIHJzSwkGUlpby8ccfM3XqVL9LERERkS6gTcHKzM4ws/VmtsHMbmql/QYzW2tmH5nZm2Y2xPtSj7wVK1ZQV1enYCUiIiJtcshgZWZRwAPAmUAWcKGZZTVb7UNginNuHPAc8BuvC/XDsmXLABSsREREpE3aMmJ1HLDBObfROVcFLAbObbyCc+4t51x56OW/gXRvy/THsmXLGD58OKmpqX6XIiIiIl1AdBvWGQRsbfS6ADjYEM7lwCutNZjZPGAeQFpaGtnZ2W2rsp1KS0s7dIx33nmH8ePHd3qdkaaj/SKdQ/0SntQv4Un9Ep7CoV/aEqyslWWu1RXNLgKmACe31u6cewR4BGDKlClu+vTpbauynbKzs2nvMQoKCigsLGTmzJnt3oe0riP9Ip1H/RKe1C/hSf0SnsKhX9oSrAqAjEav04FtzVcys1OBW4CTnXOV3pTnH11fJSIiIoerLddYfQCMNLOhZhYLzAGWNF7BzCYCDwMznXM7vS/zyFu2bBmxsbFMmDDB71JERESkizhksHLO1QDXAK8B64BnnXM5ZrbAzGaGVrsDSAL+bGarzGzJAXbXZSxbtoyJEycSFxfndykiIiLSRbTlVCDOuZeBl5st+3Gj56d6XJevamtrWb58Od/61rf8LkVERES6EM283opNmzZRXl6u04AiIiJyWBSsWrF27VoAsrKaz4MqIiIicmAKVq2oD1ZjxozxuRIRERHpShSsWpGTk0NGRga9evXyuxQRERHpQhSsWrF27VqdBhQREZHDpmDVTF1dHevWrVOwEhERkcOmYNXMli1b2L9/v4KViIiIHDYFq2b0iUARERFpLwWrZvSJQBEREWkvBatmcnJyGDBgAH369PG7FBEREeliFKya0ScCRUREpL0UrBpxzrF27VrGjh3rdykiIiLSBSlYNbJ161bKyso0YiUiIiLtomDViD4RKCIiIh2hYNWIgpWIiIh0hIJVI7m5uaSmppKSkuJ3KSIiItIFKVg1kpeXx4gRI/wuQ0RERLooBatGNm7cyLBhw/wuQ0RERLooBauQqqoq8vPzGT58uN+liIiISBelYBWyZcsW6urqNGIlIiIi7aZgFbJx40YAjViJiIhIuylYheTl5QEKViIiItJ+ClYhGzduJD4+nv79+/tdioiIiHRRClYheXl5DBs2jEBAPxIRERFpH6WIEE21ICIiIh2lYAU458jLy9P1VSIiItIhClbAzp07KSsr04iViIiIdIiCFZpqQURERLyhYIWmWhARERFvKFgRHLEyMzIzM/0uRURERLowBSuCI1aDBg0iPj7e71JERESkC1OwQlMtiIiIiDcUrEBTLYiIiIgnIj5YlZeXs337dgUrERER6bCID1abNm0C0KlAERER6bA2BSszO8PM1pvZBjO7qZX2L5nZSjOrMbNZ3pfZeTTVgoiIiHjlkMHKzKKAB4AzgSzgQjPLarZaPjAXeNrrAjtb/eSgGrESERGRjopuwzrHARuccxsBzGwxcC6wtn4F59zmUFtdJ9TYqfLy8ujVqxcpKSl+lyIiIiJdXFuC1SBga6PXBcDU9hzMzOYB8wDS0tLIzs5uz27arLS09JDH+OCDD+jXrx9vv/12p9Yin2tLv8iRp34JT+qX8KR+CU/h0C9tCVbWyjLXnoM55x4BHgGYMmWKmz59ent202bZ2dkc6hjFxcWMHz/+kOuJd9rSL3LkqV/Ck/olPKlfwlM49EtbLl4vADIavU4HtnVOOUdWbW0tmzZt0oXrIiIi4om2BKsPgJFmNtTMYoE5wJLOLevI2LZtG1VVVbpwXURERDxxyGDlnKsBrgFeA9YBzzrncsxsgZnNBDCzL5hZAXA+8LCZ5XRm0V7RVAsiIiLipbZcY4Vz7mXg5WbLftzo+QcETxF2KZpqQURERLwU0TOv5+XlERUVxeDBg/0uRURERLqBiA5WGzduZMiQIURHt2ngTkREROSgIjpY5eXl6foqERER8YyClYKViIiIeCRig1VxcTG7d+/WhesiIiLimYgNVvWfCNSIlYiIiHgl4oOVRqxERETEKxH7cbj6yUEVrERERA6surqagoICKioq/C7lkHr37s26des6tI/4+HjS09OJiYlp1/YRG6w2btxIamoqvXr18rsUERGRsFVQUEDPnj3JzMzEzPwu56BKSkro2bNnu7d3zlFUVERBQQFDhw5t1z4i9lTgJ598wogRI/wuQ0REJKxVVFSQkpIS9qHKC2ZGSkpKh0bnIjZY5ebmMmbMGL/LEBERCXuREKrqdfR7jchgVVxczI4dOxg9erTfpYiIiEg3EpHBav369QAasRIRERFPRWSwys3NBdCIlYiIiHgqIoPVunXriImJafcV/yIiInLkbN68mdGjR3PppZcybtw4Zs2aRXl5ud9ltSoip1vIzc1l5MiRREdH5LcvIiLSLtdffz2rVq3ydJ8TJkzgnnvuOeR669ev57HHHmPatGlcdtllLFy4kO9973ue1uKFiByxys3N1WlAERGRLiQjI4Np06YBcNFFF/Hee+/5XFHrIm7Iprq6mry8PP7nf/7H71JERES6lLaMLHWW5tMghOsUEBE3YpWXl0dNTY1GrERERLqQ/Px8li5dCsAzzzzDiSee6HNFrYu4YFX/iUBNtSAiItJ1jBkzhkWLFjFu3Dh2797N/Pnz/S6pVRF3KrD+5oyjRo3yuRIRERFpq0AgwEMPPeR3GYcUkSNWgwYN6tBNGkVERERaE5HBStdXiYiIdB2ZmZmsWbPG7zLaJKKClXNOwUpEREQ6TUQFqx07drBv3z4FKxEREekUERWsPvzwQwDGjh3rcyUiIiLSHUVUsMrOziYmJoapU6f6XYqIiIh0QxEXrKZOnUpCQoLfpYiIiEg3FDHBat++faxYsYIZM2b4XYqIiIj4bNeuXUydOpWJEyfy7rvverbfiJkg9L333qOuro7p06f7XYqIiIgcAbW1tURFRbXa9uabbzJ69GgWLVrk6TEjJlhlZ2cTGxvL8ccf73cpIiIiXVdrAxQXXABXXQXl5XDWWS3b584NPgoLYdaspm3Z2Yc85ObNmznjjDOYOnUqH374IUcffTRPPPFEq5f2ZGZmctlll/H6669zzTXX8IUvfIGrr76aXbt2kZCQwKOPPkpFRQU33ngj+/fvZ8KECSxdupQePXq05bs/pIg5Fajrq0RERLqu9evXM2/ePD766CN69erFwoULD7hufHw87733HnPmzGHevHncf//9rFixgjvvvJOrrrqKCRMmsGDBAmbPns2qVas8C1UQISNW9ddX3XLLLX6XIiIi0rUdbIQpIeHg7ampbRqhak1GRgbTpk0D4KKLLuK+++7je9/7Xqvrzp49G4DS0lLef/99zj///Ia2ysrKdh2/rSIiWL377ru6vkpERKQLM7ODvm4sMTERgLq6OpKTk1m1alWn1tZYRJwKrL++6otf/KLfpYiIiEg75Ofns3TpUgCeeeYZTjzxxENu06tXL4YOHcqf//xnIHhru9WrV3dqnd0+WO3cuZM//vGPnHDCCZ6eQxUREZEjZ8yYMSxatIhx48axe/du5s+f36bt/vjHP/LYY48xfvx4xo4dywsvvNCpdbbpVKCZnQHcC0QBv3fO3d6sPQ54ApgMFAGznXObvS318NXW1vLNb36T3bt3c/fdd/tdjoiIiLRTIBDgoYceOuR6mzdvbvJ66NChvPrqqy3Wmzt3LnPnzvWous8dcsTKzKKAB4AzgSzgQjPLarba5cAe59wI4G7g114X2h5PPPEEb7zxBg888AATJkzwuxxGav4nAAAgAElEQVQRERHp5tpyKvA4YINzbqNzrgpYDJzbbJ1zgfoZtp4DTrGDXVV2BLz22ms8+eSTzJ07l8suu8zPUkRERKQDMjMzWbNmTZNlX//615kwYUKTxxtvvOFThZ9ry6nAQcDWRq8LgOZ3MW5YxzlXY2Z7gRSgsPFKZjYPmAeQlpZGdjs/ctkWxcXFzJgxg9mzZ/P222932nHk8JWWlnZq30v7qF/Ck/olPEVSv/Tu3Zt9+/Yd9FN4fnjiiSdaLKutraWkpKRD+3XOUVFR0e7+bUuwau0n6dqxDs65R4BHAKZMmeI6e/qDtLQ0TbEQhrKzs9UvYUj9Ep7UL+Epkvpl06ZNVFVVkZKSEnbhqrmSkhJ69uzZ7u2dcxQVFZGcnMzEiRPbtY+2BKsCIKPR63Rg2wHWKTCzaKA3sLtdFYmIiEjYSE9Pp6CggF27dvldyiFVVFQQHx/foX3Ex8eTnp7e7u3bEqw+AEaa2VDgU2AO8I1m6ywBLgWWArOAfzrnWoxYiYiISNcSExPD0KFD/S6jTbKzs9s90uSVQwar0DVT1wCvEZxu4Q/OuRwzWwAsd84tAR4DnjSzDQRHquZ0ZtEiIiIi4ahN81g5514GXm627MeNnlcA5zffTkRERCSSdPuZ10VERESOFPPrUigz2wVs6eTDpNJsygcJC+qX8KR+CU/ql/CkfglPndkvQ5xzfQ+1km/B6kgws+XOuSl+1yFNqV/Ck/olPKlfwpP6JTyFQ7/oVKCIiIiIRxSsRERERDzS3YPVI34XIK1Sv4Qn9Ut4Ur+EJ/VLePK9X7r1NVYiIiIiR1J3H7ESEREROWIUrEREREQ80m2DlZmdYWbrzWyDmd3kdz0CZpZhZm+Z2TozyzGz//W7JvmcmUWZ2Ydm9pLftUiQmSWb2XNmlhv6d/NFv2sSMLPvhH6HrTGzZ8ysY3f9lXYxsz+Y2U4zW9No2VFm9g8z+2/oa58jXVe3DFZmFgU8AJwJZAEXmlmWv1UJUAN81zk3BjgeuFr9Elb+F1jndxHSxL3Aq8650cB41D++M7NBwHXAFOfcMQTvoav74/rjceCMZstuAt50zo0E3gy9PqK6ZbACjgM2OOc2OueqgMXAuT7XFPGcc9udcytDz0sI/pEY5G9VAmBm6cBXgd/7XYsEmVkv4EsEb3KPc67KOVfsb1USEg30MLNoIAHY5nM9Eck59w6wu9nic4FFoeeLgK8d0aLovsFqELC10esC9Ac8rJhZJjARWOZvJRJyD3AjUOd3IdJgGLAL+L/QKdrfm1mi30VFOufcp8CdQD6wHdjrnHvd36qkkTTn3HYI/s880O9IF9Bdg5W1skzzSoQJM0sC/gJc75zb53c9kc7MzgZ2OudW+F2LNBENTAIedM5NBMrw4bSGNBW6ZudcYCgwEEg0s4v8rUrCSXcNVgVARqPX6WioNiyYWQzBUPVH59zzftcjAEwDZprZZoKnzb9sZk/5W5IQ/D1W4JyrH9V9jmDQEn+dCmxyzu1yzlUDzwMn+FyTfO4zMxsAEPq680gX0F2D1QfASDMbamaxBC8sXOJzTRHPzIzg9SLrnHO/9bseCXLO3eycS3fOZRL8t/JP55z+D9xnzrkdwFYzGxVadAqw1seSJCgfON7MEkK/005BHyoIJ0uAS0PPLwVeONIFRB/pAx4JzrkaM7sGeI3gJzb+4JzL8bksCY6MXAx8bGarQst+6Jx72ceaRMLZtcAfQ/+DuBH4fz7XE/Gcc8vM7DlgJcFPOn9IGNxGJRKZ2TPAdCDVzAqAnwC3A8+a2eUEQ/D5R7wu3dJGRERExBvd9VSgiHQjZjbXzN7rwPZXmNk9B2nfbGantnf/jfZznZnd3tH9iEjXpWAlIp4ys2wz22NmcX7XAhA6jfYj4I4jcLhHgIvM7Ih/xFtEwoOClYh4JjQ/2UkEpzeZ6WsxnzsXyA3NP9SpnHMVwCvAJZ19LBEJTwpWIuKlS4B/E7zVxKWNG8zscTN7wMz+bmYlZrbMzIY3aj89dH/PvWa20MzeNrNvtXYQMxsdug/Y7tA2FxykpjOBt5ttf7GZbTGzIjO7pVlbwMxuMrO8UPuzZnZUo/ZLGm17ayunEbMJzmIvIhFIwUpEvHQJ8MfQ4ytmltas/ULgZ0AfYANwG4CZpRKcp+lmIAVYzwHmBgrNPv4P4GmCsypfCCw0s7EHqOnY0P7qt88CHiT4CdWBoeOlN1r/OoK3wTg51L6H4L1H67ddCHwTGAD0puVdHdYRvK+fiEQgBSsR8YSZnQgMAZ4NzeKeB3yj2WrPO+f+45yrIRi+JoSWnwXkOOeeD7XdB+w4wKHOBjY75/7POVcTuv/kX4BZB1g/GShp9HoW8JJz7h3nXCVwK01v5XMFcItzriDU/lNgVui+cLOAF51z74XuQ/pjWt7VoYRg4BKRCNQt57ESEV9cCrzunCsMvX46tOzuRus0DkvlQFLo+UAa3d/TOedC89K0Zggw1cwa35A4GnjyAOvvAXo2et38WGVmVtRs/381s8ZhqxZIa2Xb8mbbEjrW3gPUIiLdnIKViHSYmfUALgCizKw+PMUByWY23jm3+hC72E6j03GhGa3TD7DuVuBt59xpbSzvI+DoZsca0+hYCQRPBzbe/2XOuX8135GZbQdGNXrdo9m2hPZ9qO9XRLopnQoUES98jeCoThbB03sTCAaMd2nbJ+T+DhxrZl8LnXK7Guh/gHVfAo4OXYAeE3p8wczGHGD9lwleL1XvOeBsMzsxNBXDApr+LnwIuM3MhgCYWV8zO7fRtueY2QmhbX9Gy5u+n0zwk4EiEoEUrETEC5cC/+ecy3fO7ah/AL8DvhkKSwcUOn14PvAboIhgQFsOVLaybglwOsH7Gm4jeHrx1wRHyFrzIjDazAaGts8hGNyeJjh6tYfgDY/r3UvwfmOvm1kJwU85Tm207bUEb1a9neD1VDvr6zSzeILXiy062PcrIt2XbmkjImHHzAIEw843nXNvebC/eUCWc+76DhfXdL9JQDEw0jm3ycyuBTKcczd6eRwR6ToUrEQkLJjZV4BlwH7g+wRHlYY55/b7WlgzZnYO8CbBU4B3ERzNmuT0y1RE0KlAEQkfXyQ4RUMhcA7wtXALVSHnEjwFuQ0YCcxRqBKRehqxEhEREfGIRqxEREREPOLbPFapqakuMzOzU49RVlZGYmJipx5DDp/6JTypX8KT+iU8qV/CU2f2y4oVKwqdc30PtZ5vwSozM5Ply5d36jGys7OZPn16px5DDp/6JTypX8KT+iU8qV/CU2f2i5ltact6OhUoIiIi4hEFKxERERGPHFawMrM/mNlOM1vTaNlPzexTM1sVepzlfZkiIiIi4e9wr7F6nOAtKp5otvxu59ydnlQkIiIiR0R1dTUFBQVUVFT4XYonevfuzbp16zq0j/j4eNLT04mJiWnX9ocVrJxz75hZZruOdITt3FfBkrwqTjypjugonfEUERFprqCggJ49e5KZmYlZ8/uJdz0lJSX07Nmz3ds75ygqKqKgoIChQ4e2ax+HPUFoKFi95Jw7JvT6p8BcYB/Bm6Z+1zm35wDbzgPmAaSlpU1evHhxu4pui5Wf1XDfh5XMHx/H1AG+ffhRWlFaWkpSUpLfZUgz6pfwpH4JT92lX3r37s3w4cO7RagCqK2tJSoqqkP7cM6Rl5fH3r17myyfMWPGCufclENt70WwSiN4CwoH/BwY4Jy77FD7mTJliuvM6Rbq6hxf/MUr9D+qF3+7elq3edN0B/qYcnhSv4Qn9Ut46i79sm7dOsaMGeN3GZ7p6IhVvdZ+LmbWpmDV4XNkzrnPnHO1zrk64FHguI7u0wuBgPGVzBhWF+xl+ZZWB9BEREREPNXhYGVmAxq9/Dqw5kDrHmnTBkXTJyGGR9/Z6HcpIiIichgyMzMpLCxssXzJkiXcfvvtPlTUNod18ZGZPQNMB1LNrAD4CTDdzCYQPBW4GbjC4xrbLS7KuOj4IfzurQ1sKixjaKpuPyAiItKVzZw5k5kzZ/pdxgEd1oiVc+5C59wA51yMcy7dOfeYc+5i59yxzrlxzrmZzrntnVVse1z8xSHEBAI8/q9NfpciIiIizWzevJnRo0dz6aWXMm7cOGbNmkV5eTkA999/P5MmTeLYY48lNzcXgMcff5xrrrnGz5IPqtt/XK5fz3imjUjhP5t1nZWIiMiB/OzFHNZu2+fpPrMG9uIn54w95Hrr16/nscceY9q0aVx22WUsXLgQgNTUVFauXMnChQu58847+f3vf+9pfZ0hIiZ4GpqaxJaiMg73E5AiIiLS+TIyMpg2bRoAF110Ee+99x4A5513HgCTJ09m8+bNfpV3WLr9iBVAZmoC5VW17CqppF+veL/LERERCTttGVnqLM2nRKp/HRcXB0BUVBQ1NTVHvK72iIgRqyEpwYvWNxeV+1yJiIiINJefn8/SpUsBeOaZZzjxxBN9rqj9IiJYZaYkALC5qMznSkRERKS5MWPGsGjRIsaNG8fu3buZP3++3yW1W0ScChyU3IPogLFFwUpERCTsBAIBHnrooSbLGl9TNWXKFLKzswGYO3cuc+fOPXLFHaaIGLGKjgqQcVSCTgWKiIhIp4qIYAUwJCWBzYUasRIREQknmZmZrFkTNjdt6bCICVaZKYlsKSrXlAsiIiLSaSImWA1JSaC0soaisiq/SxEREZFuKmKCVWboPoG6gF1EREQ6S+QEq/q5rAp1AbuIiIh0jogJVoOSexClKRdERESkE0VMsIqNDjAouQebNOWCiIhI2MvMzKSwsLDF8iVLlnD77bd7cowLL7yQcePGcffdd3uyP4iQCULrDUlJ0IiViIhIFzZz5kxmzpzZpnVramqIjm496uzYsYP333+fLVu2eFle5IxYAQxNTWRTYZmmXBAREWnN9OktHwsXBtvKy1tvf/zxYHthYcu2Nti8eTOjR4/m0ksvZdy4ccyaNYvy8uDZpfvvv59JkyZx7LHHkpubC8Djjz/ONddcc8D9zZ07lxtuuIEZM2bwgx/8gLKyMi677DK+8IUvMHHiRF544QUATj/9dHbu3MmECRN499132/oTOqSIClZDUhIpqaihuLza71JEREQkZP369cybN4+PPvqIXr16sTAU5lJTU1m5ciXz58/nzjvvbPP+PvnkE9544w3uuusubrvtNr785S/zwQcf8NZbb/H973+fsrIylixZwvDhw1m1ahUnnXSSZ99LRJ0KrL8Z86aiMvokxvpcjYiISJgJ3Y+vVQkJB29PTT14+0FkZGQwbdo0AC666CLuu+8+AM477zwAJk+ezPPPP9/m/Z1//vlERUUB8Prrr7NkyZKGYFZRUUF+fj49evRoV62HElHBalCf4A9xx94KnysRERGRembW6uu4uDgAoqKiqKmpafP+EhMTG5475/jLX/7CqFGjmqzT+CbPXoqoU4EpicEOKiqt9LkSERERqZefn8/SpUsBeOaZZzjxxBM92/dXvvIV7r///obrqz/88EPP9t2aiApWfRJiMIPCUt3WRkREJFyMGTOGRYsWMW7cOHbv3s38+fM92/ett95KdXU148aN45hjjuHWW2/1bN+tiahTgdFRAfokxFJUphErERGRcBEIBHjooYeaLGt8qm7KlClkh67fmjt3LnPnzj3gvh6v/5RiSI8ePXj44YdbrJeZmcmaNWvaW/IBRdSIFUBqUiyFJRqxEhEREe9F1IgVBK+z0oiViIhIeGjvyNFtt93Gn//85ybLZs6cyYIFC7wqrV0iL1glxbJ22z6/yxAREQkLzrkWn8rrCm655RZuueWWJstKSko6vN+OTiIegacC4yjUpwJFRESIj4+nqKhIdyQJcc5RVFREfHx8u/cReSNWibHsq6ihqqaO2OiIy5UiIiIN0tPTKSgoYNeuXX6X4omKiooOhSIIhs309PR2bx9xwSq1Z2guq7JKBvTunFlXRUREuoKYmBiGDh3qdxmeyc7OZuLEib7WcNhDNmb2BzPbaWZrGi07ysz+YWb/DX3t422Z3kkJ3cqmSHNZiYiIiMfacy7sceCMZstuAt50zo0E3gy9DkspScERK11nJSIiIl477GDlnHsH2N1s8bnAotDzRcDXOlhXp0lN0oiViIiIdA5rzycBzCwTeMk5d0zodbFzLrlR+x7nXIvTgWY2D5gHkJaWNnnx4sXtLLttSktLSUpKarJsf41j/hvlzB4Vy5lDYzr1+NK61vpF/Kd+CU/ql/CkfglPndkvM2bMWOGcm3Ko9Y7oxevOuUeARwCmTJnipk+f3qnHy87OpvkxnHPEv/0qvdPSmT59TKceX1rXWr+I/9Qv4Un9Ep7UL+EpHPrFq/kGPjOzAQChrzs92q/nzIyURM1lJSIiIt7zKlgtAS4NPb8UeMGj/XaK1KRYXWMlIiIinmvPdAvPAEuBUWZWYGaXA7cDp5nZf4HTQq/DVkqS7hcoIiIi3jvsa6yccxceoOmUDtZyxKQkxrJuu+4XKCIiIt6KyHu6pPaMo6i0SvdGEhEREU9FZLBKSYylqraOfRU1fpciIiIi3UhEBqvU0OzrRfpkoIiIiHgoIoNVSv3s62X6ZKCIiIh4JzKDVaJGrERERMR7ERmsUnsGR6wKNZeViIiIeCgig9VRCfXBSiNWIiIi4p2IDFbRUQH6JMRo9nURERHxVEQGK9Ds6yIiIuK9yA1WibG6xkpEREQ8FbHBKjUpTtdYiYiIiKciNlj1SYyhuLza7zJERESkG4ncYJUQS3F5FXV1ul+giIiIeCOig1Wdg30VGrUSERERb0RusEqMAWCPTgeKiIiIRyI3WIUmCd2t+wWKiIiIRyI+WBWXK1iJiIiINyI2WB2VqBErERER8VbEBqvkhOA1VppyQURERLwSscEqKS6a6ICxW6cCRURExCMRG6zMjD6JsbrGSkRERDwTscEKoE9CjK6xEhEREc9EeLCK1TxWIiIi4pmID1Y6FSgiIiJeiexglRjL7jKNWImIiIg3IjtYJcRQXF6Fc7oRs4iIiHRcRAeroxJjqalzlFTW+F2KiIiIdAMRHayS629ro9OBIiIi4oFoL3dmZpuBEqAWqHHOTfFy/17rE5p9fXd5FYNTEnyuRkRERLo6T4NVyAznXGEn7NdzfUL3C9yjTwaKiIiIByL6VGCf0KnAPZokVERERDxgXn4izsw2AXsABzzsnHukWfs8YB5AWlra5MWLF3t27NaUlpaSlJR0wPayasfVb5Zz4ehYvpIZ06m1yOcO1S/iD/VLeFK/hCf1S3jqzH6ZMWPGirZc4uT1qcBpzrltZtYP+IeZ5Trn3qlvDAWtRwCmTJnipk+f7vHhm8rOzuZgx6ircwT++TIpAwYzffqoTq1FPneofhF/qF/Ck/olPKlfwlM49IunpwKdc9tCX3cCfwWO83L/XgsELHRbG50KFBERkY7zLFiZWaKZ9ax/DpwOrPFq/50lOSFGwUpEREQ84eWpwDTgr2ZWv9+nnXOverj/TnFUYix7NI+ViIiIeMCzYOWc2wiM92p/R0pyQixbd5f7XYaIiIh0AxE93QLAUbrGSkRERDwS8cEqOTGGPWXVuhGziIiIdFjEB6s+CbFU1dZRXlXrdykiIiLSxUV8sDoqNPv6bs2+LiIiIh0U8cEqOXQj5uJyfTJQREREOibig9VRoRsx79YF7CIiItJBER+sknUjZhEREfFIxAerlNCIVZGClYiIiHRQxAer3j1iiA4YRaWVfpciIiIiXVzEB6tAwEhJiqVQwUpEREQ6KOKDFUBKYhyFpToVKCIiIh2jYAWk9ozTiJWIiIh0mIIVkJoUS5FGrERERKSDFKyAvklx7Cqt1P0CRUREpEMUrICUpFiqauooqazxuxQRERHpwhSsgNSkOAAKS3SdlYiIiLSfghWfBytNEioiIiIdoWCFRqxERETEGwpWBD8VCGjKBREREekQBSvgqMRYzGCXplwQERGRDlCwAqKjAvRJiNX9AkVERKRDFKxCUnW/QBEREekgBasQ3S9QREREOkrBKkT3CxQREZGOUrAK0f0CRUREpKMUrEJSk+IorayhorrW71JERESki1KwCqmfy2qXJgkVERGRdlKwCmmYfV3XWYmIiEg7eRaszOwMM1tvZhvM7Cav9nukNNwvUNdZiYiISDtFe7ETM4sCHgBOAwqAD8xsiXNurRf7PxJSe2rEqrNUVNdSXF5NaWUN+6tqWb+7FtbvZH9VLeVVteyvrqWiupbqWkdNbR01dY7aOkd1XR21tY6aOkdNXR21dY6aWocL7dfqv9rnx7LQ0vplwa9GwCA6YAQCRpQZUY2e13+NjjICZkQFCH0NrWd2wG2jAhAVCDTdpn6fAWtYFrBguxkNywMGZqH1Q20NNdW3BT5/Hmi0rYX2FxV6bo1/CCIi4htPghVwHLDBObcRwMwWA+cCvgarCddfD8nJTRdecAFcdRWUl8NZZzUsHuAcizftprTuIjjuB1BYCLNmtdzp/PkwezZs3QoXX9yy/bvfhXPOgfXr4YorWrb/6Edw6qmwahVcf33L9l/+Ek44Ad5/H374w5bt99wDEybAG2/AL37Rsv3hh2HUKHjxRbjrrpbtTz4JGRnwpz/Bgw+2bH/uOUhNhccfDz6ae/llXI8elN1zP/bss1TV1lFVU0dVbR3VtY7f/OAhivdXccrfn2Tyx/+iptZR54JRqCI6jrkX/AyAa+++gWlbVtOn0a739OjF/K8Hv+cfvLOISZ/mBkMDweCwq3dfFsy+mYAZ1734ACO3b2hSWn5KOrd//QYcjpv/ejcZhQUNbQ5Y338Yt59+JbXOcftff0PavkKccw3rrBw0ml+fPBeAB//6S/rs39dk//8aMp77p10IwOPP/oT4mqYh/M3hx/Ho1PMAWPz0TdQCtUB1qP2l0Sfx1KSvEl9dweN//mnLH/2xp/LcsafSp3wvD/7tVy3an5p4Fi+N+RID9u3i7pda9u1jU8/j7aOPZ9juAha8fP/nyZPg04dPupBlwyYxakceP3jtIZq2wqPTL+Y771QxPj+Ha958vFkr3HXWfP47cATHbVjJZdl/bGisb//1177D1r4ZTFu3lAvfe67F9r+YcxM7k/txyupszv33khb1/+Tin7I3sTdnLH+VM5a/1qL9B5f9isrYeL72/gtM/yi7Rfv1V94NwOy3/8QX1/27SVtVTBw3Xn47AJe88SSTNqxs0r4voRc/viT43vz2K48ydkvTX127evfltguD781rlvyOEdvymrQXpKZz56zvAvC95+4ivdF7D2DDwOH8buY1ANzyzC/pu3dXk/acIVk8eua3AVjwxE/oVf75e69nXR2PHT2FRacEf9/85rGbiGv23ls6+ngWnzwbgHsf/k6Ln81b46bzty+eS1xVBb/5v5tbtL8y+Su8OuUMepftZcFTP23R/sLxM/nn+Bn0K97JLX9q+d7800nn837WCWTsyud7z9/dov2JL1/EipGTGbFtA9e++ECL9ke/cjlrMo/hmM1r+PZrj7Vov/+cq9kwcAST/7uCS/75VIv2O8/7Dlv7DuaEte8z+90/t2j/xezge+/Lq9/ia/9+sUX7rRf9hL2JvTlz+aucueL1Fu3f/3+/DL73lr7Alz96G4CetXV8GBU86XPdFb8FYM7bz3JCbtP3XmV0LN8PvfcufeNJJud92KR9b0Ivbr34pwBc8crvGZvf/L2Xys/nBN971y55gJHbm773tqamc8f/3BCs8y+/bfJ7D+C/A4Zz/8yrg9/n4l/Sd29hk/acwVk8fOa3APj5kz+ld3nT33srhk9k0anB994dj91EXE3Tszrvjz6exSdfAMB9D99Ac/8cd3LDe++O/2v5N+2VyafzSui99/Onftai/W/Hn9Pw3vvRn25v0b74pFm8n3UCV548nG9OHdKi3Q9eBatBwNZGrwuAqc1XMrN5wDyAtLQ0srOzPTp8646traW4uLjJsp2ffMK27GwCFRWMa9YWMNi6YyfZ2dnE7N3L2GbtAJ/m5LArO5u4nTsZ00r71o8/pqhnT3rk5zOqlfYtq1ezJzqapA0bGNFK+8aVK9lXVUWvNWsY1kr7huXLKS0ups/q1QxppX39smXs376dlI8/JqOV9nVLl1KZl0ffnBwGtdKe869/Ud27N/1zc+m7p5jKWkdVLcFHneN/b3+VrbVxXPBBLmd/urfJtlEB46Mtn5EUY8TUVdEjyhEVDVEWIBCAutgA102Mo666gpO3xJBeHEUAQiMyUNMrit+fnkDAYPjWGHqXNz1TndHX+NUXowAYsdxIKms6StM/zbj9hGD70e8bCdXN2gdFMWhGcGRyzNIo4kK/FOujVdrQGIaflkCdgwn/iiJ2XxT1ucsBSZkxZHypB3UOjv5HgKiqKByO0H+cOjiavlPjqXOQ8WLTfQPMyIimz8Q4AhWOga80a3dw0qBoeh4TS499saQlBhqW16/zxYHRJIyKoVdRNH17BBqW13+d3C+K2MFR9I0NkBzXtA1g7FEQ6A/p1dArtvFPJjgKOCrZEZ/qGLYPkqIbbxk0omcdsb1rGZxYS0Ir7YMTaohNqqF/fC09our+f3t3Hh5Xdd9//P3VrNplS7ZsS17kHbPZRmACAWRwgBAgS8P6QCFtf6aENaHJL4H+CEnbtGnJRpImoYEmD0lxUkoWKIRdEJYYvLDYyAZjvMibbHmTLGsZzfn9MWNZy8jyckd3pPm8nkePZ+73Ll/pjO2vzj33HHrvMTYSIxxtZ0Swg4jF+xw/JtRGfriNkkAsZXxsuI32sFEcTB0fF04UG0UpjjeLd8ULU8SjOd3iOZ194nk5nV3xglTxwMF4XqBvvKDb8Xkpji/sFo/mxHvEXY6jMBCjItIGBtFAnHC85/FFoRgV0cTxkZy+P5uSYCIezmlLGR+RPL4w1l+8g4poGyMj7SnjpeFEvDzckTJeloyPDqc+viySiEv1WcsAACAASURBVJdFUh8/OtzO/mgbZf2cvzzcQTzaRmk/8bGRGJFoB6WhzpTxikgHRdEORvYTr4x20B4OMCJ4MO5wWPL1+Gji16eSFMdbwHXFi1PE83LiXfHC4EDxeJ94QeBgvCDQN14YPBjPy0kV7zxkvDh0MB4NuD6fvZJu8ZSfneT5wzmp22Zk8vjCWOp4aTI+MhJLGR8VTsS3r3+f2v0f0tzcnPbaYiDW/Tf2oz6J2WXABc65v0m+vxY4zTl3S3/HVFdXuyVLlhzztQ+ltraWmpqaw95//r21HD+uiB9ePTd9SWWoLXv2s3zDblZu3sOqLU3UbdnL5j2tXfFwIIcJpXlMKs1jYmk+40fkMqY4SnlRlDHFUUYVRAgGDm/I3pG2iwwOtUtmUrtkJrVLZkpnu5jZUudc9UD7edVjVQ+M7/a+Etjs0bkHTWl+dqwXGOuM8/amPSxdt4vlG3exbP1utu5NFFHBHGPKqAJOrRrJzDFFzBxbyLTRBYwtziWQo3E8IiIih+JVYfUGMM3MqoBNwJXA1R6de9CUFURYs73Z7zQ855xj9bYmXlnTyKtrdrD4w500t8UAGD8yl9OqRjJnQglzJ4xg5thCIsGAzxmLiIgMTZ4UVs65mJndDDwFBIAHnXMrvTj3YBpVGOG1tY1+p+GJ9licVz/YwdPvbuOZd7d1TXw6qTSPS2eP48wpZZxaNYLRhVGfMxURERk+vOqxwjn3BPCEV+fzQ8WIXPbs76CptYPCaMjvdI5YW6yT5+saeHLFVl5Y1UBTW4y8cID5M0ZzzoxRnDm1jIqSXL/TFBERGbY8K6yGgwkj8wDYsLOF48cV+5zN4Xt/WxOL3tjIo8vq2dXSwcj8MB8/cQwXHD+GM6eWEQ3p1p6IiMhgUGHVzYHCauMQKKxaOzp57K3NLHpjI0vX7yIUMM6fNYYrTh3PGVNKD/sJPREREfGOCqtuxnfrscpUTa0d/GrxBh54+UO2N7UxeVQ+d110HJ+eW9G1LI+IiIj4Q4VVN8W5IUryQhlZWLV2dPKLV9fx4xc/YHdLB2dNK+N7V8zmjCmlWs5EREQkQ6iw6mXCyDw27NzvdxpdOjrj/PeSeu577n227m3lnOmjuOP86ZxUWTLwwSIiIjKoVFj1Mn5kHit7LdXih3jc8fg7W/jO06tZ19jCKRNH8P0rZzNvcqnfqYmIiEg/VFj1MmFkHk+t2Epn3Pky07hzjmfrGvj206tZtbWJmWMKeeC6as6dOVq3/ERERDKcCqteJozMIxZ3bNmzn8oReYN2Xeccr6xp5N+eXs1bG3dTVZbP96+czSUnjSNHS8mIiIgMCSqseuk+l9VgFFZtsU6eXrmNh15bz+vrdjKuOMq3/uJE/mJupaZMEBERGWJUWPXSfS4rpqTnGvvaYrz6QSMvrG7gjyu2snNfO5UjcrnnkllcNW+C1uoTEREZolRY9TK2OEowxzyfcqF+VwtPrdxG7eoGFq/dSXtnnPxwgHNmjOKKUydw1tQy3fITEREZ4lRY9RIM5FAxIpf1jcdeWHXGHY+9tZnfLNnIqx8kFneeOrqA686YyPwZo6meNJJwULf7REREhgsVVilMGJmXuBV4DNbt2MeXHnmLN9btYsLIPL74sel8cvY4Jpbme5SliIiIZBoVVimMH5nHk+9sOerjH1laz9//7h1CgRy+fdnJfGZuhaZKEBERyQIqrFKYMDKPXS0d7G3toCgaOqJj/7hiK1965C1Oryrlu1fMZkxxNE1ZioiISKbRAJ8UJnZ/MvAILNuwi9sWLefkyhIevP5UFVUiIiJZRoVVCuOPorCq39XC3/xiCWOKozxwXTW5YU2ZICIikm10KzCFCaWJwmrtjn2Htb9zjrt/v5LWjk7+58YzKC2IpDM9ERERyVDqsUqhKBpi6ugC/rx252Ht/9TKrTy/qoEvLJhOVZme+hMREclWKqz6cfa0USxe20hrR+ch92tui3HPH97luLFFfO7MSYOTnIiIiGQkFVb9OGt6GW2xOK9/eOheq+8+8x7bmlr5p0+foLX9REREspwqgX6cXlVKOJjDS+9t73efD7Y384tX13HlqeOZO2HEIGYnIiIimUiFVT9ywwFOmzSSP72/o999vvXkKiLBHL74sRmDmJmIiIhkKhVWh3DWtDJWb2ti657WPrHFaxt5+t1t3FgzhVGFegpQREREVFgd0tnTRwHw0vs9bwfG445vPlHH2OIof/3RyX6kJiIiIhlIhdUhzBxTyKjCSJ/bgQ/9eT1v1e/h786foYlARUREpIsKq0MwM86aVsZL721nTUMzAM+v2sbXH1vJeTNH8+k5FT5nKCIiIplEhdUArvvIJMzgovv+xD8/WcfN/7Wc48cV84Or55CTY36nJyIiIhnkmAsrM7vHzDaZ2ZvJr4u8SCxTnDy+hKe/cDbnTB/FT19cy4i8MA9cV01eWKsBiYiISE9eVQffdc7d69G5Ms7owij3X3sKL763namjCxhdFPU7JREREclA6nY5TGZGzYzRfqchIiIiGcycc8d2ArN7gOuBvcAS4A7n3K5+9l0ILAQoLy8/ZdGiRcd07YE0NzdTUFCQ1mvIkVO7ZCa1S2ZSu2QmtUtmSme7zJ8/f6lzrnqg/Q6rsDKzZ4ExKUJ3AX8GdgAO+AdgrHPurwY6Z3V1tVuyZMmA1z4WtbW11NTUpPUacuTULplJ7ZKZ1C6ZSe2SmdLZLmZ2WIXVYd0KdM4tOMyL/gfw+OHsKyIiIjLcePFU4Nhubz8NrDjWc4qIiIgMRV4MXv9XM5tN4lbgOuAGD84pIiIiMuQc8+D1o76w2XZgfZovU0Zi/JdkFrVLZlK7ZCa1S2ZSu2SmdLbLROfcqIF28q2wGgxmtuRwBprJ4FK7ZCa1S2ZSu2QmtUtmyoR20ZI2IiIiIh5RYSUiIiLikeFeWN3vdwKSktolM6ldMpPaJTOpXTKT7+0yrMdYiYiIiAym4d5jJSIiIjJoVFiJiIiIeGTYFlZmdqGZrTazNWb2Fb/zETCz8Wb2gpnVmdlKM7vN75zkIDMLmNlyM9OyVBnCzErM7BEzW5X8e/MRv3MSMLMvJP8NW2FmD5tZ1O+cspGZPWhmDWa2otu2kWb2jJm9n/xzxGDnNSwLKzMLAD8CPg7MAq4ys1n+ZiVADLjDOXcccDpwk9olo9wG1PmdhPTwfeCPzrmZwMmofXxnZhXArUC1c+4EIABc6W9WWevnwIW9tn0FeM45Nw14Lvl+UA3Lwgo4DVjjnFvrnGsHFgGf9DmnrOec2+KcW5Z83UTiP4kKf7MSADOrBD4B/MzvXCTBzIqAs4EHAJxz7c653f5mJUlBINfMgkAesNnnfLKSc+4lYGevzZ8EfpF8/QvgU4OaFMO3sKoANnZ7X4/+A88oZjYJmAMs9jcTSfoe8GUg7nci0mUysB34z+Qt2p+ZWb7fSWU759wm4F5gA7AF2OOce9rfrKSbcufcFkj8Mg+MHuwEhmthZSm2aV6JDGFmBcD/ALc75/b6nU+2M7OLgQbn3FK/c5EegsBc4MfOuTnAPny4rSE9JcfsfBKoAsYB+WZ2jb9ZSSYZroVVPTC+2/tK1FWbEcwsRKKo+pVz7lG/8xEAzgQuNbN1JG6bn2tmv/Q3JSHx71i9c+5Ar+4jJAot8dcC4EPn3HbnXAfwKHCGzznJQdvMbCxA8s+GwU5guBZWbwDTzKzKzMIkBhb+weecsp6ZGYnxInXOue/4nY8kOOe+6pyrdM5NIvF35XnnnH4D95lzbiuw0cxmJDedB7zrY0qSsAE43czykv+mnYceKsgkfwCuS76+Dvj9YCcQHOwLDgbnXMzMbgaeIvHExoPOuZU+pyWJnpFrgXfM7M3ktjudc0/4mJNIJrsF+FXyF8S1wOd8zifrOecWm9kjwDISTzovJwOWUclGZvYwUAOUmVk98DXgX4DfmNlfkyiCLxv0vLSkjYiIiIg3huutQBGRw2Jmr5jZnCPY/ztm9rfpzElEhi4VViIyKMxsnZntN7Pmbl/jfM7pEqDJObc8+f4eM3Nmdlm3fYLJbZOSm/4NuCt5e05EpAcVViIymC5xzhV0+zqip3WTEzJ66W+Bh3pt2wl8I7mCQx/JuXFWAZd6nIuIDAMqrETEd2Z2aXLttd1mVmtmx3WLrTOz/2tmbwP7kj1I483sUTPbbmaNZvbDbvv/VXJdvV1m9pSZTeznmmHgXODFXqE/Au3AoZ6MrCUxU72ISA8qrETEV2Y2HXgYuB0YBTwBPNbrVttVJAqZEhKT/T4OrAcmkVhVYVHyXJ8C7gQ+kzzXn5LnTmUaEHfO1ffa7oD/B3wtOe9aKnUk1u4TEelBhZWIDKbfJXuldpvZ75LbrgD+1zn3THLCxXuBXHpOunifc26jc24/ibVAxwFfcs7tc861OudeTu53A/DPzrk651wM+CYwu59eqxKgKVWSzrk/kFhO5m/6+T6akseLiPSgwkpEBtOnnHMlya8Di6OOI9H7BIBzLk5irc/u63t2X/tzPLA+WTj1NhH4/oHijcR4KSP1WqG7gMJD5Pr3wF1ANEWsENCCyCLShworEfHbZhIFEdA1Q/94YFO3fbpPuLcRmNDPQPaNwA3dircS51yuc+7VFPu+n7xcygXanXPPAGuAz6cIHwe8dahvSkSykworEfHbb4BPmNl5yTFNdwBtQKpiCOB1YAvwL2aWb2ZRMzszGfsJ8FUzOx7AzIq7T53QXfK247PAOYfI7S7gyym2nwM8OcD3JSJZSIWViPjKObeaxBN4PwB2AJeQmJahvZ/9O5P7TCWxZEU9iXFaOOd+C3wLWGRme4EVwMcPcfmfklhmqb/cXiFRyHVJLuw6C/hdyoNEJKtpSRsRyWpm9jJwy4FJQg9j/28DHzjn/j29mYnIUKTCSkRERMQjuhUoIiIi4hFPCyszC5jZcjN73MvzioiIiAwFXvdY3UZiRmIRERGRrOPZgqZmVkliyYl/Ar440P5lZWVu0qRJXl0+pX379pGfn5/Wa8iRU7tkJrVLZlK7ZCa1S2ZKZ7ssXbp0h3Nu1ED7eblS/PdIzPfS70zGZrYQWAhQXl7Ovffe6+Hl+2pubqagoCCt15Ajp3bJTGqXzKR2yUxql8yUznaZP3/++oH38qiwMrOLgQbn3FIzq+lvP+fc/cD9ANXV1a6mpt9dPVFbW0u6ryFHTu2SmdQumUntkpnULpkpE9rFqzFWZwKXmtk6EqvMn2tmv/To3CIiIiJDgieFlXPuq865SufcJOBK4Hnn3DVenFtERERkqPByjJWIiIgMMx0dHdTX19Pa2up3KgMqLi6mru7YJieIRqNUVlYSCoWO6njPCyvnXC1Q6/V5j8VL720nHMzh9MmlfqciIiIypNTX11NYWMikSZMwM7/TOaSmpiYKC/t9hm5AzjkaGxupr6+nqqrqqM4x7Gdeb2mPccvDy/nO0+/5nYqIiMiQ09raSmlpacYXVV4wM0pLS4+pd27YF1a/W76ZPfs72L2/3e9UREREhqRsKKoOONbvdVgXVs45fv7qhwDsaunwORsREREZ7oZ1YVW3M85725oZVxxlT0sHzjm/UxIREZFhbFgXVk+v66A0P8wVp06gvTPO/o5Ov1MSERGRYWzYTrewvnEfb23v5JZzqygvigCwu6WDvPCw/ZZFRETS6uuPreTdzXs9PeescUV87ZLjD7nPunXruPjii1mxYgUA9957L83Nzdxzzz2e5uKFYdtjta6xhRFR45rTJ1KSFwZgV4sGsIuIiEj6DNvum3Omj+Lec3IZXRSlJC8xydceDWAXERE5agP1LMkw7rECyEk+MnmgsNq9X4WViIjIUBMMBonH413vM3kW+GFdWB1Qkpu4FbhbPVYiIiJDTnl5OQ0NDTQ2NtLW1sbjjz/ud0r9Gra3Ars72GOlMVYiIiJDTSgU4u6772bevHlUVVUxc+ZMv1PqV1YUVtFQgGgoRz1WIiIiQ9Stt97Krbfe6ncaA8qKW4GQuB24W08FioiISBplT2GVF1KPlYiIiKRV1hRWxbkhPRUoIiIiaZU1hVVJXkjzWImIiEhaZU1hNSIvrJnXRUREJK2yprAqzkvcCnTO+Z2KiIiIDFOeFVZmFjWz183sLTNbaWZf9+rcXijJDdMei9PaER94ZxEREZGj4GWPVRtwrnPuZGA2cKGZne7h+Y+JJgkVERGRA9ra2liwYAGzZ8/m17/+tWfn9WyCUJe4x9acfBtKfmXMfbeS3GRh1dLB2OJcn7MREREZompq+m67/HL4/OehpQUuuqhv/PrrE187dsBnP9szVlvrfY5JsViMYDB1qbN8+XI6Ojp48803Pb2mpzOvm1kAWApMBX7knFvcK74QWAiJdX9q0/jDBGhubu66xrrGTgBqX32DbaWBtF5XDq17u0jmULtkJrVLZsqmdikuLqapqanrfW5nZ599Yq2tdDQ1QUtLynhHayuxpiasuZlor/j+bufuz/r167n88stZvDhRVtx33300Nzdz55139tivs7OTs846i3nz5vHnP/+Ziy66iKuuuorbb7+djRs3AvCtb32LKVOmcPXVV9PY2MhJJ53EQw89xOTJk7vO09raetTt62lh5ZzrBGabWQnwWzM7wTm3olv8fuB+gOrqaleTqur1UG1tLQeuUb5lL996409Mmj6LmhPHpvW6cmjd20Uyh9olM6ldMlM2tUtdXR2FhYUHN/zpT332CQJRgMLCfuP0Ey/ss3dfBQUF5OTkdOURiUTo6OjomRfQ1NREIBCgpaWFl19+GYCrr76aL33pS3z0ox9lw4YNXHDBBdTV1fHggw9y7733plzQORqNMmfOnMPIrK+0rBXonNttZrXAhcCKAXYfFAfHWGkuKxERkeHsiiuu6Hr97LPP8u6773a937t3b48eOK95VliZ2SigI1lU5QILgG95df5jVZIbBtCyNiIiIkNMMBgkHj/4VH9ra+sh98/Pz+96HY/Hee2118jNHZzx1V4+FTgWeMHM3gbeAJ5xzvXtX/NJNJRDOJijpwJFRESGmPLychoaGmhsbKStrS3l7bv+nH/++fzwhz/seu/1YPXevHwq8G3g6G5IDgIzY0ReiN371GMlIiIylIRCIe6++27mzZtHVVUVM2fOPOxj77vvPm666SZOOukkYrEYZ599Nj/5yU/SlmtaxlhlqpLcsHqsREREhqBbb72VW2+9dcD9ej/NV1ZWlnKeqpqamrQ8gJA1S9pAclkbjbESERGRNMmyHqsQG3a2+J2GiIiIHKObbrqJV155pce2G264gRtvvNGnjBKyq7DKC/F2vXqsREREjoRzDjPzO40efvSjH/XZ5sU0ComFZI5eVt0KHJEXZleLxliJiIgcrmg0SmNj4zEXHEOBc47Gxkai0ehRnyOreqyK80K0xeK0dnQSDWlZGxERkYFUVlZSX1/P9u3b/U5lQK2trcdUFEGikKysrDzq47OqsOo+SeiYYhVWIiIiAwmFQlRVVfmdxmGpra096qVovJJVtwIPLmuj24EiIiLivewsrDTlgoiIiKRBVhVWRdFEYbVHCzGLiIhIGmRVYVUQSQwp29cW8zkTERERGY6yqrDKV2ElIiIiaZRVhVVhNFFYNamwEhERkTTIqsIqEswhkGPqsRIREZG0yKrCyszIDwfY19bpdyoiIiIyDGVVYQVQGA3R1KoeKxEREfFe1hVW+ZGAbgWKiIhIWnhSWJnZeDN7wczqzGylmd3mxXnToSASpFmFlYiIiKSBV2sFxoA7nHPLzKwQWGpmzzjn3vXo/J7JjwR1K1BERETSwpMeK+fcFufcsuTrJqAOqPDi3F4riAR1K1BERETSwpxz3p7QbBLwEnCCc25vr9hCYCFAeXn5KYsWLfL02r01NzdTUFDQY9sD77SxsrGT79TkpfXa0r9U7SL+U7tkJrVLZlK7ZKZ0tsv8+fOXOueqB9rPq1uBAJhZAfA/wO29iyoA59z9wP0A1dXVrqamxsvL91FbW0vva9TuXcmbjfV9tsvgSdUu4j+1S2ZSu2QmtUtmyoR28eypQDMLkSiqfuWce9Sr83rtwK1Ar3vqRERERLx6KtCAB4A659x3vDhnuhREg8Qd7O/QJKEiIiLiLa96rM4ErgXONbM3k18XeXRuTx1YiFlTLoiIiIjXPBlj5Zx7GTAvzpVuhQcKq9YYowt9TkZERESGlSyceT1RWGm9QBEREfFaFhZWAQCa2jp8zkRERESGm6wrrAojIUA9ViIiIuK9rCusDvRYafZ1ERER8VrWFVYFyTFWTSqsRERExGPZV1hFDwxeV2ElIiIi3sq6wio3FCDHVFiJiIiI97KusDIz8iNBmlpVWImIiIi3sq6wgoPrBYqIiIh4KSsLq/xIUEvaiIiIiOeysrAqUGElIiIiaZC1hZVuBYqIiIjXsrKwyo8E1GMlIiIinsvKwqogEtKSNiIiIuK5LC2sAjS1ahFmERER8VZ2FlbRIPvaO3HO+Z2KiIiIDCNZWVjlR4J0xh1tsbjfqYiIiMgwkpWFVddCzJp9XURERDzkWWFlZg+aWYOZrfDqnOlyoLDSlAsiIiLiJS97rH4OXOjh+dImP1lYacoFERER8ZJnhZVz7iVgp1fnS6cCFVYiIiKSBublk3FmNgl43Dl3Qj/xhcBCgPLy8lMWLVrk2bVTaW5upqCgoM/2tXs6+cZrrdw+N8Ls0cG05iB99dcu4i+1S2ZSu2QmtUtmSme7zJ8/f6lzrnqg/Qa1qnDO3Q/cD1BdXe1qamrSer3a2lpSXaOyoRlee5Gq6cdRM7sirTlIX/21i/hL7ZKZ1C6ZSe2SmTKhXbLyqcDCqG4FioiIiPeysrDqGryu6RZERETEQ15Ot/Aw8Boww8zqzeyvvTq31/JCAUDTLYiIiIi3PBtj5Zy7yqtzpVtOjlEQCdKshZhFRETEQ1l5KxAgPxKguU0LMYuIiIh3sriwCrJPPVYiIiLioawtrAojQZo0xkpEREQ8lLWFVaLHSoWViIiIeCdrC6uCSFDTLYiIiIinsrawKoyGNEGoiIiIeCprC6ui3CB79+upQBEREfFO9hZW0RBNbTE6494tQi0iIiLZLXsLq9wQAE2t6rUSERERb2RvYZVciHnvfo2zEhEREW9kb2GV7LHaqx4rERER8Uj2FlbRZGGlAewiIiLikewtrHKTtwLVYyUiIiIeyd7CqqvHSmOsRERExBtZW1gV52mMlYiIiHgrawurgnAQM42xEhEREe9kbWGVk2MURoLs1XqBIiIi4pGsLawgMeWCeqxERETEK54VVmZ2oZmtNrM1ZvYVr86bTkXRkMZYiYiIiGeCXpzEzALAj4CPAfXAG2b2B+fcu16cP10SCzHrVuDR6uiM09rRSWvHgT8Tr/d3dBLrjNMRd8Q648TijlinIxaPE+t0rKjvYNPi9cQ6HR2dcTrjjlhyzUbnHC65fKODbq8PxKFrdcdksPt+B5gl/+y2wbrFDrxLvO51jHUddch9u28/uL91Xbd7Dl3bu+eVcl/rk3v3a9Ejh17fU7dzpMrrULkbsHJrjNYVWxJnSJV7f9caaN8UPz/6/bn21y6H2Lefn99wsW5PJys27fE7DelF7ZJZRhdFGF0Y9TsNwKPCCjgNWOOcWwtgZouATwK+Flazb78dSkp6brz8cvj856Glha9/92baOuKwqPhg/PrrE187dsBnP9v3pDfeCFdcARs3wrXX9o3fcQdccgmsXg033NA3/vd/DwsWwJtvwu23941/85twxhnw6qtw551949/7HsyeDc8+C//4j33jP/0pzJgBjz0G3/523/hDD8H48fDrX8OPf0ync3TE4nR0Oto74zz/Dz+iIVLIxD/8hllP/zZRDCULo8644+Zr/oE9FuYzf/49F9X9qc/pr7z6XwD4P4sf5bwPXifSLdYajHD95V8H4Jaf3sWZ69/qceyu3CJu/HTie/7yiz9n7qZVPeJbCsv4wiV/B8Ddz97PrIa1PeLrSiu46+O34pzjn578AVU7N/WIvzt6Mt9YsBCA7z52L2ObdvSIL6uYyb+ecz0AP/7tNxmxf2+P+CsTT+YHZ14FwM9/8zWisbYe8eemnMZ/zPsMAIv+q2+n7eMzz+KXcz9BtKOVn//3PX3ij5y4gEdOXMCIlj38+Hf/3Cf+yzkX8fhxZzN273a++3jftv2P0z7Nc1PnMbmxnm8+9cM+8R+ccSWvTJrNrG1rufu5+/vEXz37On74Zhtz6+v48ku/6BP/xnkLebd8Mmeue5NbXl3UJ37nBTeztrSS89Ys5v+8/ts+8S9cfAdbikZxcd1LXLP8iR4xB9z4qa+yK6+Yz77zLJ9959k+x19/2T20hqJcs+x/uXjVwJ+97np89l552PPP3tqRFdx54S0AfPOPP2Cyh5+9mcBTw/yz969nX8eyyuN8+ezB0X32ZgLNyfhw/ezB0Pl37ysfn8nfnjOlT9wPXhVWFcDGbu/rgXm9dzKzhcBCgPLycmpraz26fGondnaye/fuHtsa3nuPzbW15LS2MirWQXuMHvtsXbWKrbW1hPbs4fhexwJsWrmS7bW1RBoaOC5FfOM779BYWEjuhg3MSBFf/9Zb7AoGKVizhqkp4muXLWNveztFK1YwOUV8zZIlNO/ezYi33mJiivjqxYvZv2ULBcvfZnzj7kRRFIeYg1jc8f2f1fJBXhmnLl3OpR82Eu/V0/PNJ+rYlVfMFe9vZ8y+NnIMAgY5BpEAHD8ijkXh+NIAZbk5WDJmJP68bW6EUI5xypYgE7bndPUmGEZnOId/PTuX1pYWTt0cYtTeQCKevHasOMDPzs8DYOqGEMX7Az1ymzQqyH9ekIeZMXVNiIJYz/iEiiAPJI+f/k6QPHrGKyeFqLogDwfMWhoksr1nvLwqxPTzE/GTXg0Q2tMzXjolzAkLEvHjnw+Q09YzXjwtxOzzEtef9uTB2IEfviXWNgAADXNJREFU8eUzwpx2bh6B/TlMfabnsQBXzQxxZk0u4T3tTKntG7/muDDnnJ1L7vYok1/uFXdw/awwH5uXS2F9lMl/DtCrafmrE8JcNDuXER9EqFrSN7+/nO741OwoZXVhJr2V0zMI/O1JYRonRxmbF2Liir6jCG48OcKeiijjc8JMWNU3ftPsCPvKokzuCDF+Td/4LXOitBZFmdEUYvy6g/EDKdw6J0IsEuHEHUEqN/U9/ra5iTL+lE1BKrbldOseg1g4pys+78MgFTt7Hl9acDBevSrI2L0948UjA13x2W8HGLW/Zzy/7GD8xCUBSjp6xsOjD8ZnvhKgoNcojJzyYFd86gs5RIMH4/HOOB8ZFyQnGZ/4VA7B9p7Hn1UZJC8Zr3is789m/vggpXMjBNscFX/sG//YxCAVcyNE90aoeK5v/ONVIabMjVCwI0zFi33jl04JccLcCCM2h6l4tW/8M1NDVJ8UYdS6EBWv941fPiPEWTMijM0PUbG8b/zq40JsnxRhQjBExTt94395fJhd4yJMdiEq6vrG/+qEMM1lEaa3hqhI8dlbeFKE1qIIs/YGqViX+rMbi0Q4aXuQiuRnL94ZJyeQeN3js9fQ8/ih/NkDhs5nr3k9tbUbaW5uTnttMRBzve+hHM1JzC4DLnDO/U3y/bXAac65W/o7prq62i1ZsuSYr30otbW11NTU9Bv/xmPv8us3NrDyGxemNY902LO/gzUNTdTv2s+m3fvZvHs/m3e3smlX4nVTW99bnHnhAKMLI4wqTHSZjiqMdHWfHthemh+mKDdENNT3P3avDNQu4g+1S2ZSu2QmtUtmSme7mNlS51z1QPt51WNVD4zv9r4S2OzRudOmKDfIvvbEeKBg4NjG8cfjjmUbdpEXDlJRkktRbrDHWJ1jPXfd1r28/uFOXv9wJ2/X72HT7v099inJCzGuOJcJpXl8ZEop40qijCnOZXRhJPFVFKUg4lVzi4iISCpe/U/7BjDNzKqATcCVwNUenTttDixr09QaY0R++KjP0xbr5Iu/fov/fWdL17YZ5YV85eMzqZkx6qgKrM6445U1O3hyxVaeq9tGQ1PivnbliFzmTCjh6nkTOG5sIRNG5jG2OJd8FU0iIiK+8+R/Y+dczMxuBp4CAsCDzrmVXpw7nYpyDy5rc7SFVVNrBzc8tJRXP2jkix+bztTRBWzc2cLDr2/gcz9/g49OLeOrF83k+HHFA58MqN/Vwn8vqee/l2xk855W8sMBzpkxivNmlid7onKPKk8RERFJP8+6OZxzTwB9H7fIYEXRxLd/LFMufP5Xy3j9w51894qT+fScyq7tnzuzil8tXs/3n3ufi3/wMn8xt5LbF0yjckRen3Psbmnnxfe288jSel5ek3hi46NTy7jrE7M477jRaR3vJCIiIt7J6vtH3XusjsaHO/bxp/d38KULZvQoqgDCwRw+d2YVn5lbyb+/sIb/fGUdjyytZ8qofOZNLiUcyGFvawdrt+/j7frdxB1UlORy67nTuKy6MmUBJiIiIpktqwur4gOF1VEua/PosnpyDD57SmW/+xTnhvjqRcdxzekTeXLFFl79oJHH3twMlhjjNaY4yi3nTuPs6aOYPb6EQM4wm91QREQki2R1YXUsPVbxuOPRZZs4c2oZ5UUDz/Y6fmQeC8+ewsKzM2MCMxEREfFedi/CfAxjrBZ/uJNNu/cfsrdKREREsktWF1b54SA5dnQ9Vo8uq6cgEuT8WWPSkJmIiIgMRVldWOXkGIXR0BGPsWppj/HEO1u46MQx5Ib1xJ6IiIgkZHVhBYnZ1/e2HtmtwOdXNbCvvbPPk4AiIiKS3VRYHUWP1eK1O8kPBzh10og0ZSUiIiJDkQqraOiIx1gtWb+LORNGHPP6giIiIjK8ZH1lUJQbPKKnAptaO1i9dS+nTFRvlYiIiPSkwioaYs8R3ApcviExS3q1bgOKiIhILyqsco/sVuCS9bvIMZg9viSNWYmIiMhQpMIqGqKlvZOOzvhh7b90/U5mjCmiMBpKc2YiIiIy1Kiwyk3Mvt50GFMuxDrjLN+wm2qNrxIREZEUVFhFD38h5lVbm2hp79T4KhEREUlJhdURLMS8dP0uAD0RKCIiIillfWFVfKCwOowpF5as38WYoigVJbnpTktERESGoGMurMzsMjNbaWZxM6v2IqnBdGCM1WH1WK3bySmTRmBm6U5LREREhiAveqxWAJ8BXvLgXIOuJDcMQOO+9kPut21vK5v3tDJ3gm4DioiISGrBYz2Bc64OGLK9OKMKI4QDOdTvajnkfss37AY0f5WIiIj0z5xz3pzIrBb4O+fckkPssxBYCFBeXn7KokWLPLl2f5qbmykoKBhwv6+81EJlYQ43z4n2u89vVrfz1LoOfrwgj3BgaBaRmeJw20UGl9olM6ldMpPaJTOls13mz5+/1Dk34JCnw+qxMrNngTEpQnc5535/uEk55+4H7georq52NTU1h3voUamtreVwrjHzw9fZ3tRGTc1Z/e7zk/de4/iKTs4/76MeZpidDrddZHCpXTKT2iUzqV0yUya0y2EVVs65BelOxE8TRuaxdN0unHMpb2l2xh3v1O/hL06p9CE7ERERGSqyfroFSBRWTW2xfhdjfr+hiX3tnRpfJSIiIofkxXQLnzazeuAjwP+a2VPHntbgGj8yD4ANO1MPYH9TA9dFRETkMBxzYeWc+61zrtI5F3HOlTvnLvAiscE0sTRRWK1v7Kew2rib4twQVWX5g5mWiIiIDDG6FQiMHzFAj9XG3Zw8vmTITikhIiIig0OFFZAfCVJWEGZjisJqX1uM97Y16TagiIiIDEiFVdKEkXkpe6zert9D3MEcFVYiIiIyABVWSRNG5qUcY7Vswy4ATlZhJSIiIgNQYZU0YWQeW/bspz0W77H9xfe2M2tsESPzwz5lJiIiIkOFCquk8SPziDvYvHt/17a9rR0sXb+LmhmjfMxMREREhgoVVkkTSxNTKazvNs7qlfd30Bl31MwY7VdaIiIiMoSosEqakGKS0BdWN1AYDTJ3gsZXiYiIyMBUWCWNLowQDuZ0TbngnOPF97Zz9rRRBAP6MYmIiMjAVDEk5eRYYsqF5JOBdVua2La3jXM0vkpEREQOkwqrbiaMzGNd4z6cc9S+1wBAzXQVViIiInJ4gn4nkEmOH1fE86sauPSHr9DcFmPW2CJGF0X9TktERESGCPVYdXPT/Kn8w6dOoC3WyYc79rFgVrnfKYmIiMgQoh6rbqKhANeePpFr5k1g1dYmqsry/U5JREREhhAVVimYGceNLfI7DRERERlidCtQRERExCMqrEREREQ8osJKRERExCMqrEREREQ8osJKRERExCPmnPPnwmbbgfVpvkwZsCPN15Ajp3bJTGqXzKR2yUxql8yUznaZ6JwbcDkW3wqrwWBmS5xz1X7nIT2pXTKT2iUzqV0yk9olM2VCu+hWoIiIiIhHVFiJiIiIeGS4F1b3+52ApKR2yUxql8ykdslMapfM5Hu7DOsxViIiIiKDabj3WImIiIgMGhVWIiIiIh4ZtoWVmV1oZqvNbI2ZfcXvfATMbLyZvWBmdWa20sxu8zsnOcjMAma23Mwe9zsXSTCzEjN7xMxWJf/efMTvnATM7AvJf8NWmNnDZhb1O6dsZGYPmlmDma3otm2kmT1jZu8n/xwx2HkNy8LKzALAj4CPA7OAq8xslr9ZCRAD7nDOHQecDtykdskotwF1fichPXwf+KNzbiZwMmof35lZBXArUO2cOwEIAFf6m1XW+jlwYa9tXwGec85NA55Lvh9Uw7KwAk4D1jjn1jrn2oFFwCd9zinrOee2OOeWJV83kfhPosLfrATAzCqBTwA/8zsXSTCzIuBs4AEA51y7c263v1lJUhDINbMgkAds9jmfrOScewnY2WvzJ4FfJF//AvjUoCbF8C2sKoCN3d7Xo//AM4qZTQLmAIv9zUSSvgd8GYj7nYh0mQxsB/4zeYv2Z2aW73dS2c45twm4F9gAbAH2OOee9jcr6abcObcFEr/MA6MHO4HhWlhZim2aVyJDmFkB8D/A7c65vX7nk+3M7GKgwTm31O9cpIcgMBf4sXNuDrAPH25rSE/JMTufBKqAcUC+mV3jb1aSSYZrYVUPjO/2vhJ11WYEMwuRKKp+5Zx71O98BIAzgUvNbB2J2+bnmtkv/U1JSPw7Vu+cO9Cr+wiJQkv8tQD40Dm33TnXATwKnOFzTnLQNjMbC5D8s2GwExiuhdUbwDQzqzKzMImBhX/wOaesZ2ZGYrxInXPuO37nIwnOua865yqdc5NI/F153jmn38B95pzbCmw0sxnJTecB7/qYkiRsAE43s7zkv2nnoYcKMskfgOuSr68Dfj/YCQQH+4KDwTkXM7ObgadIPLHxoHNupc9pSaJn5FrgHTN7M7ntTufcEz7mJJLJbgF+lfwFcS3wOZ/zyXrOucVm9giwjMSTzsvJgGVUspGZPQzUAGVmVg98DfgX4Ddm9tckiuDLBj0vLWkjIiIi4o3heitQREREZNCpsBIRERHxiAorEREREY+osBIRERHxiAorEREREY+osBIRERHxiAorEREREY/8f+H93810ySHcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot results\n",
    "\n",
    "fig,axes = plt.subplots(3,1, figsize=(10,10))\n",
    "axes[0].plot(tsim, xsim[:,0], \"k\", label='p')\n",
    "axes[0].plot(tsim, xref[0]*np.ones(np.shape(tsim)), \"r--\", label=\"p_ref\")\n",
    "axes[0].set_title(\"Position (m)\")\n",
    "\n",
    "axes[1].plot(tsim, xsim[:,2]*360/2/np.pi, label=\"phi\")\n",
    "axes[1].plot(tsim, xref[2]*360/2/np.pi*np.ones(np.shape(tsim)), \"r--\", label=\"phi_ref\")\n",
    "axes[1].set_title(\"Angle (deg)\")\n",
    "\n",
    "axes[2].plot(tsim, usim[:,0], label=\"u\")\n",
    "axes[2].plot(tsim, uref*np.ones(np.shape(tsim)), \"r--\", label=\"u_ref\")\n",
    "axes[2].set_title(\"Force (N)\")\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid(True)\n",
    "    ax.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
